home *** CD-ROM | disk | FTP | other *** search
/ BBS in a Box 2 / BBS in a box - Trilogy II.iso / Files / Education / M / MathPad 2.15 / examples / Navigation < prev    next >
Encoding:
Text File  |  1993-10-08  |  2.8 KB  |  95 lines  |  [TEXT/MPAD]
  1. -- calculate great circle routes
  2.  
  3. places = {Anchorage, BuenosAires, CapeTown, Honolulu, London,
  4.        Moscow, NewYork, Singapore, Sydney, Tokyo, Wellington}
  5.  
  6. here=London; there=Anchorage
  7.  
  8. az    := heading(here,there):;            az:-15.6
  9. alpha := angle(here,there):; distance(alpha):4457.0
  10.  
  11. -- For a nice background, use the world map picture from the system 7 distribution Scrapbook File. Paste it into the PLOT and comment out the following line.
  12. plot place(places)
  13.  
  14. plot place({here})
  15.  
  16. s := rotation(here[1],here[2],90):
  17. inc := multiply(rotation(incangle,0,az),rotation(0,0,-az)):
  18.  
  19. plot {step,c[1]}; plot place({there})
  20.  
  21. step = r:={1,0,0} when X=0,
  22.        r:=transform(inc,r), b:=transform(s,r), c:=geog(b),
  23.        c[2]
  24.  
  25. place(L)[i] = {L[i,2],L[i,1]}  -- swap for plot. X=lon Y=lat
  26.  
  27. increment=300; incangle=increment/LenPerDeg
  28. Xmin=-180;  Xmax=180;  Xsteps=trunc(alpha/incangle)
  29. Ymin= -89;  Ymax= 89
  30.  
  31. -- locations are given in {+N latitude, +E longitude} degrees.
  32. Anchorage  ={ 61,-150}
  33. BuenosAires={-35,- 58}
  34. CapeTown   ={-34,  18}
  35. Honolulu   ={ 21,-157}
  36. London     ={ 52,   0}
  37. Moscow     ={ 56,  38}
  38. NewYork    ={ 41,- 74}
  39. Singapore  ={  1, 103}
  40. Sydney     ={-34, 151}
  41. Tokyo      ={ 36, 140}
  42. Wellington ={-41, 174}
  43.  
  44. -- the great circle heading from A to B
  45. -- is the azimuth of B-A at A
  46. heading(A,B) = azimuth(head(A,B))
  47.  
  48. -- the arc length of an angle A
  49. distance(A) = A*LenPerDeg
  50.  
  51. LenPerDeg = radius*π/180
  52.  
  53. -- the radius of earth (pick distance units)
  54. radius = 3960
  55.  
  56. -- azimuth is measured east of north
  57. -- in the spherical system V[2] is south, V[3] is east
  58. azimuth(V) = atan2(V[3],-V[2])
  59.  
  60. -- B-A at A
  61. head(A,B) = rotate(A[2],-A[1],-90,cart(B)-cart(A))
  62.  
  63. -- angle between locations A and B
  64. angle(A,B) = acos(sum(cart(A)[i]*cart(B)[i],i,1,3))
  65.  
  66. -- cartesian coordinates for unit vector at A={lat,lon}
  67. cart(A) = {sin(90-A[1])*cos(A[2]),sin(90-A[1])*sin(A[2]),cos(90-A[1])}
  68.  
  69. -- {lat,lon} of cartesian unit vector A
  70. geog(A) = {90-acos(A[3]),atan2(A[2],A[1])}
  71.  
  72. atan2(y,x) = 0 when x=0 and y=0,
  73.              atan(y/x) when x>=0,
  74.              atan(y/x)+180 when y>=0,
  75.              atan(y/x)-180
  76.  
  77. -- the vector V in a system rotated by yaw, pitch, and roll
  78. rotate(yaw,pitch,roll,V) = transform(rotation(yaw,pitch,roll),V)
  79.  
  80. rotation(y,p,r) = {{cos(p)*       cos(y),
  81.                     cos(p)*       sin(y),
  82.                    -sin(p)                            },
  83.                    {sin(p)*sin(r)*cos(y)-cos(r)*sin(y),
  84.                     sin(p)*sin(r)*sin(y)+cos(r)*cos(y),
  85.                     cos(p)*sin(r)                     },
  86.                    {sin(p)*cos(r)*cos(y)+sin(r)*sin(y),
  87.                     sin(p)*cos(r)*sin(y)-sin(r)*cos(y),
  88.                     cos(p)*cos(r)                     }}
  89.  
  90. transform(A,V)[i] = sum(A[i,k]*V[k],k,1,count(V)) dim [count(A)]
  91.  
  92. multiply(A,B)[i,j] = sum(A[i,k]*B[k,j],k,1,count(B)) dim [count(A),count(B[1])]
  93.  
  94. transpose(A)[i,j] = A[j,i] dim [count(A[1]),count(A)]
  95.